function printout_fields_of_eigenmodes(m, n, CONSTS, span, plot_data)

    a = CONSTS.a;
    
    p_m_vec = func_p_m(m, CONSTS, span, false);
    coef_m_mat = func_coef_m_of_eigenmodes(m, CONSTS, p_m_vec);

%     [p_m_vec, coef_m_mat] = pmvec_coefmmat_of_eigenmodes_analytically(m, CONSTS, span); 
%     [p_m_vec, coef_m_mat] = pmvec_coefmmat_of_eigenmodes(m, CONSTS, span);
%     fprintf('p_mn  %f \n', p_m_vec);
    
    p_n = p_m_vec(n);
    coef_n = (coef_m_mat(:,n));
    rho_high = 10;
    num_rho_points = 5000; %200
    rho_step = rho_high/num_rho_points;
    rho_1_vec = (0:rho_step:a)';
    rho_2_vec = (a:rho_step:rho_high)';
    r1_vec = rho_1_vec./a;
    r2_vec = rho_2_vec./a;

    [Hrho_in, Hphi_in, Hz_in, Erho_in, Ephi_in, Ez_in] = func_HE_inside(r1_vec, m, CONSTS, coef_n, p_n);
    [Hrho_out, Hphi_out, Hz_out, Erho_out, Ephi_out, Ez_out] = func_HE_outside(r2_vec, m, CONSTS, coef_n, p_n);

    if(plot_data)
        
        clear = @(x)clear_small_values(x, 1e-5);
    
        % "-" due to H
        figure; plot(r1_vec, clear(real(Erho_in)),  'b.-', ...
                     r1_vec, clear(imag(Ephi_in)),     'r.-', ...
                     r1_vec, clear(imag(Ez_in)),       'm.-', ...
                     r2_vec, clear(real(Erho_out)), 'b.-', ...
                     r2_vec, clear(imag(Ephi_out)),    'r.-',...
                     r2_vec, clear(imag(Ez_out)),      'm.-', ...
                     [1, 1], real([Erho_in(size(Erho_in,1)), Erho_out(1)]), 'b:', ...
                     [0, rho_high/a], [0, 0], 'k--');
        title('field E of eigenmode'); xlabel('\rho/a'); ylabel('E');
        legend('Re(E_\rho)', 'Im(E_\phi)', 'Im(E_z)', 'location', 'SouthEast');

        % +
        figure; plot(r1_vec, clear(imag(Erho_in)),  'b.-', ...
                     r1_vec, clear(real(Ephi_in)),     'r.-', ...
                     r1_vec, clear(real(Ez_in)),       'm.-', ...
                     r2_vec, clear(imag(Erho_out)), 'b.-', ...
                     r2_vec, clear(real(Ephi_out)),    'r.-',...
                     r2_vec, clear(real(Ez_out)),      'm.-', ...
                     [1, 1], imag([Erho_in(size(Erho_in,1)), Erho_out(1)]), 'b:', ...
                     [0, rho_high/a], [0, 0], 'k--');
        title('field E of eigenmode'); xlabel('\rho/a'); ylabel('E');
        legend('Im(E_\rho)', 'Re(E_\phi)', 'Re(E_z)', 'location', 'SouthEast');

        % Hphi, Hrho - discontinuous
        figure; plot(r1_vec, clear(imag(Hrho_in)), 'b.-', ...
                     r1_vec, clear(real(Hphi_in)), 'r.-', ...
                     r1_vec, clear(real(Hz_in)),   'm.-', ...
                     r2_vec, clear(imag(Hrho_out)), 'b.-', ...
                     r2_vec, clear(real(Hphi_out)), 'r.-', ...
                     r2_vec, clear(real(Hz_out)),   'm.-', ...
                     [0, rho_high/a], [0, 0], 'k--');
        title('field H of eigenmode'); xlabel('\rho/a'); ylabel('H');
        legend('Im(H_\rho)', 'Re(H_\phi)', 'Re(H_z)', 'location', 'SouthEast');

        % +
        figure; plot(r1_vec, clear(real(Hrho_in)), 'b.-', ...
                     r1_vec, clear(imag(Hphi_in)), 'r.-', ...
                     r1_vec, clear(imag(Hz_in)),   'm.-', ...
                     r2_vec, clear(real(Hrho_out)), 'b.-', ...
                     r2_vec, clear(imag(Hphi_out)), 'r.-', ...
                     r2_vec, clear(imag(Hz_out)),   'm.-', ...
                     [0, rho_high/a], [0, 0], 'k--');
        title('field H of eigenmode'); xlabel('\rho/a'); ylabel('H');
        legend('Re(H_\rho)', 'Im(H_\phi)', 'Im(H_z)', 'location', 'SouthEast');

    end
    
    
    
end

